MSA Cluster Modeling
Robert J. DiNapoli
9/16/2021
Load Packages
library(here)
library(mclust)
library(raster)
library(rgdal)
library(sp)
library(RColorBrewer)
library(PCAtools)
library(sf)
library(tmap)
library(terra)
library(rgeos)
library(spData)
library(knitr)Import Data
#import raw data
dat_raw <- read.csv(here("data/CommunityData-raw-2015.csv"))
#select columns of interest, removed co2_hhs
dat <- dat_raw[,c(4:12,14:21)]
#convert to dataframe
dat <- as.data.frame(unclass(dat))
#set rownames to city names
rownames(dat)=dat_raw$ME
# remove any records that have NAs
dat = na.omit(dat)
# convert to z-scores
dat <- as.matrix(scale(dat, center = TRUE, scale = TRUE))Run Cluster GMM
Assumptions: - MSAs form clusters characterized by a multivariate distribution - Model forms: shape, volume, orientation - GMM fits a series of models with different forms and numbers of clusters - Models with highest probability and fewest parameters selected as most optimal - Based on Bayesian Information Criterion (BIC)
#run mclust on 1 to 20 clusters
dat_mc <- Mclust(dat, G=c(1:20))
#print summary of model fit
summary(dat_mc)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVE (ellipsoidal, equal orientation) model with 10 components:
##
## log-likelihood n df BIC ICL
## -16206.34 886 485 -35704.23 -35902.49
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9 10
## 50 59 53 62 150 129 51 66 142 124Model comparison
#plot BIC scores of different models
plot(dat_mc, what='BIC')
#zoom in on best-fitting models
plot(dat_mc, what='BIC', ylim=c(-37000,-35000))
#model type VVE chosen as best, with similar BIC scores for 9-12 clusters
#to do additional model comparisons, generate and compare VVE models with 9-12 clusters
dat_mc9 <- Mclust(dat, G=9, modelNames="VVE")
dat_mc10 <- Mclust(dat, G=10, modelNames="VVE")
dat_mc11 <- Mclust(dat, G=11, modelNames="VVE")
dat_mc12 <- Mclust(dat, G=12, modelNames="VVE")
#extract BIC scores
BICs<-c(dat_mc9$BIC[1],dat_mc10$BIC[1],dat_mc11$BIC[1],dat_mc12$BIC[1])
#plot BIC scores
plot(c(9:12),BICs, pch=16, ylab="BIC Score",xlab="Number of Clusters", type='o')
#calculate change in BIC score, since in this method the goal to to maximize BIC
#we calculate change from highest scoring model
delta_BIC <- max(BICs) - BICs
#calculate BIC weights
w_BIC <- round(exp(-0.5*delta_BIC)/sum(exp(-0.5*delta_BIC)), digits=3)
#make table of results
results_table <- cbind.data.frame(clusters=c(9:12),BIC=BICs,delta_BIC,weight=w_BIC)
#order by delta
results_table <- results_table[order(delta_BIC),]
rownames(results_table)<-NULL
#print table
kable(results_table[1:4], caption = "Model Comparison")| clusters | BIC | delta_BIC | weight |
|---|---|---|---|
| 10 | -35704.23 | 0.00000 | 1 |
| 12 | -35731.46 | 27.23133 | 0 |
| 11 | -35811.06 | 106.82569 | 0 |
| 9 | -35818.70 | 114.46901 | 0 |
Extract and map cluster classifications
#extract classifications, e.g., which MSA belongs to which cluster, write to CSV
write.csv(dat_mc10$classification,file=here("results/cluster_assignment.csv"))
data<-read.csv("results/cluster_assignment.csv")
#import shapefile
msa_Boundary <-readOGR(here("data/MSA"),"tl_2015_us_cbsa") ## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\dinap\Desktop\Research\sustainable_communities\data\MSA", layer: "tl_2015_us_cbsa"
## with 929 features
## It has 12 fields
## Integer64 fields read as strings: ALAND AWATER#merge them
merged <- merge(msa_Boundary,data,by.x="NAME",by.y="X")
#remove any cases MSAs with no cluster assignments
merged_clean <- merged[!is.na(merged$x),]
#convert x to factor
merged_clean$Cluster <- as.factor(merged_clean$x)
#convert to sf object
cluster_sf <- st_as_sf(merged_clean)
#make us map
us_map <- tm_shape(us_states) + tm_borders()
#combine them
cluster_map <- tm_shape(cluster_sf) + tm_fill(col="Cluster", palette="Spectral")
#plot them
us_map + cluster_map
Examine how clusters score in different metrics
#new data frame
dat_clust <- as.data.frame(dat)
#add names
dat_clust$ME <- rownames(dat)
rownames(dat_clust) <- NULL
#merge in cluster assignments
dat_clust <- merge(dat_clust,data,by.x="ME",by.y="X")
#calculate cluster means for different variables
c1m <- data.frame(one=colMeans(subset(dat_clust, x=="1")[2:18]))
c2m <- data.frame(two=colMeans(subset(dat_clust, x=="2")[2:18]))
c3m <- data.frame(three=colMeans(subset(dat_clust, x=="3")[2:18]))
c4m <- data.frame(four=colMeans(subset(dat_clust, x=="4")[2:18]))
c5m <- data.frame(five=colMeans(subset(dat_clust, x=="5")[2:18]))
c6m <- data.frame(six=colMeans(subset(dat_clust, x=="6")[2:18]))
c7m <- data.frame(seven=colMeans(subset(dat_clust, x=="7")[2:18]))
c8m <- data.frame(eight=colMeans(subset(dat_clust, x=="8")[2:18]))
c9m <- data.frame(nine=colMeans(subset(dat_clust, x=="9")[2:18]))
c10m <- data.frame(ten=colMeans(subset(dat_clust, x=="10")[2:18]))
cluster_means <- cbind.data.frame(c1m,c2m,c3m,c4m,c5m,c6m,c7m,c8m, c9m, c10m)
cluster_means$variable <- row.names(cluster_means)
rownames(cluster_means) <- NULL#boxplots of scores for clusters
par(mfrow=c(2,5), mar=c(2.5,5,2.5,3))
boxplot(subset(dat_clust, x=="1")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 1")
boxplot(subset(dat_clust, x=="2")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 2")
boxplot(subset(dat_clust, x=="3")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 3")
boxplot(subset(dat_clust, x=="4")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 4")
boxplot(subset(dat_clust, x=="5")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 5")
boxplot(subset(dat_clust, x=="6")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 6")
boxplot(subset(dat_clust, x=="7")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 7")
boxplot(subset(dat_clust, x=="8")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 8")
boxplot(subset(dat_clust, x=="9")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 9")
boxplot(subset(dat_clust, x=="10")[2:18], horizontal=T, las=1, ylim=c(-4,4), main="Cluster 10")
par(mfrow=c(1,1))#dotplots of cluster mean values
par(mfrow=c(2,5))
dotchart(cluster_means$one, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 1', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$two, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 2', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$three, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 3', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$four, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 4', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$five, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 5', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$six, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 6', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$seven, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 7', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$eight, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 8', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$nine, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 9', pch=16)
abline(v=0, lwd=2)
dotchart(cluster_means$ten, labels=cluster_means$variable,xlim=c(-2,2), main='Cluster 10', pch=16)
abline(v=0, lwd=2)
par(mfrow=c(1,1))Which MSAs best characterize their cluster?
| Cluster | Name |
|---|---|
| 1 | Goldsboro, NC |
| 2 | Sierra Vista-Douglas, AZ |
| 3 | Lexington-Fayette, KY |
| 4 | Tampa-St. Petersburg-Clearwater, FL |
| 5 | Texarka, TX-AR |
| 6 | Zanesville, OH |
| 7 | Louisville/Jefferson County, KY-IN |
| 8 | Watertown, SD |
| 9 | Pama City, FL |
| 10 | Williamsport, PA |
Which cluster scores the best on the different metrics?
| Metric | Cluster |
|---|---|
| AQI_Good | eight |
| Bachelor_Over_25 | three |
| Per_Poverty | eight |
| Gini | eight |
| Per_Sev_Hous | eight |
| Xstreamlengthimpaired | two |
| Per_Avg_Land_Cov | ten |
| poor_health_percent | eight |
| Z_Water_Index | six |
| Index_Black | five |
| Index_Asian | four |
| Index_Latino | three |
| GHG_Percap | four |
| UNEMPLOY | eight |
| FOOD_INS | eight |
| VIO_CRIME | eight |
Examine cluster patterns with PCA
Note: adjusting values for variables so that high=‘good’ and low=‘bad’ for all variables
dat_clust_adj <- dat_clust
dat_clust_adj$Per_Poverty <- dat_clust$Per_Poverty * -1
dat_clust_adj$Gini <- dat_clust$Gini * -1
dat_clust_adj$Per_Sev_Hous <- dat_clust$Per_Sev_Hous * -1
dat_clust_adj$Xstreamlengthimpaired <- dat_clust$Xstreamlengthimpaired * -1
dat_clust_adj$poor_health_percent <- dat_clust$poor_health_percent * -1
dat_clust_adj$Z_Water_Index <- dat_clust$Z_Water_Index * -1
dat_clust_adj$Index_Black <- dat_clust$Index_Black * -1
dat_clust_adj$Index_Asian <- dat_clust$Index_Asian * -1
dat_clust_adj$Index_Latino <- dat_clust$Index_Latino * -1
dat_clust_adj$GHG_Percap <- dat_clust$GHG_Percap * -1
dat_clust_adj$UNEMPLOY <- dat_clust$UNEMPLOY * -1
dat_clust_adj$FOOD_INS <- dat_clust$FOOD_INS * -1
dat_clust_adj$VIO_CRIME <- dat_clust$VIO_CRIME * -1
cluster1a <- data.frame(subset(dat_clust_adj, x=="1"))
cluster2a <- data.frame(subset(dat_clust_adj, x=="2"))
cluster3a <- data.frame(subset(dat_clust_adj, x=="3"))
cluster4a <- data.frame(subset(dat_clust_adj, x=="4"))
cluster5a <- data.frame(subset(dat_clust_adj, x=="5"))
cluster6a <- data.frame(subset(dat_clust_adj, x=="6"))
cluster7a <- data.frame(subset(dat_clust_adj, x=="7"))
cluster8a <- data.frame(subset(dat_clust_adj, x=="8"))
cluster9a <- data.frame(subset(dat_clust_adj, x=="9"))
cluster10a <- data.frame(subset(dat_clust_adj, x=="10"))Perform PCA on individual clusters
#PCA cluster 1
c1_PCA <- pca(cluster1a[2:18],transposed = T) #rotate table and only use numeric columns
#scree plot
screeplot(c1_PCA, title="Scree plot Cluster 1") ## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
#elbow method of estimating important components
findElbowPoint(c1_PCA$variance) ## PC5
## 5#biplot for PC1 and PC2
biplot(c1_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 1')## Warning: Removed 1 rows containing missing values (geom_segment).## Warning: Removed 1 rows containing missing values (geom_label_repel).
#loadings for components
plotloadings(c1_PCA, title="Cluster 1")## -- variables retained:## Per_Poverty, GHG_Percap, Z_Water_Index, Index_Asian, FOOD_INS, Per_Sev_Hous, Per_Avg_Land_Cov, UNEMPLOY
#PCA cluster 2
c2_PCA <- pca(cluster2a[2:18],transposed = T)
screeplot(c2_PCA)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
findElbowPoint(c2_PCA$variance)## PC6
## 6biplot(c2_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 2')## Warning: Removed 2 rows containing missing values (geom_segment).## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c2_PCA, title="Cluster 2")## -- variables retained:## Per_Avg_Land_Cov, Per_Poverty, UNEMPLOY, Index_Asian, non_migration, Z_Water_Index, AQI_Good, Gini, Index_Latino, FOOD_INS
#PCA cluster 3
c3_PCA <- pca(cluster3a[2:18],transposed = T)
screeplot(c3_PCA)
findElbowPoint(c3_PCA$variance)## PC4
## 4biplot(c3_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 3')## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c3_PCA, title="Cluster 3")## -- variables retained:## Per_Poverty, non_migration, Index_Black, Index_Asian, Bachelor_Over_25, Per_Avg_Land_Cov, Per_Sev_Hous, AQI_Good, Index_Latino, GHG_Percap
#PCA cluster 4
c4_PCA <- pca(cluster4a[2:18],transposed = T)
screeplot(c4_PCA)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
findElbowPoint(c4_PCA$variance)## PC4
## 4biplot(c4_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 4')## Warning: Removed 2 rows containing missing values (geom_segment).## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c4_PCA, title="Cluster 4")## -- variables retained:## Xstreamlengthimpaired, AQI_Good, Index_Latino, Bachelor_Over_25, Per_Sev_Hous, Index_Black, VIO_CRIME, Per_Avg_Land_Cov, Gini
#PCA cluster 5
c5_PCA <- pca(cluster5a[2:18],transposed = T)
screeplot(c5_PCA)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
findElbowPoint(c5_PCA$variance)## PC6
## 6biplot(c5_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 5')## Warning: Removed 3 rows containing missing values (geom_segment).## Warning: Removed 3 rows containing missing values (geom_label_repel).
plotloadings(c5_PCA, title="Cluster 5")## -- variables retained:## FOOD_INS, VIO_CRIME, non_migration, Xstreamlengthimpaired, Index_Black, GHG_Percap, Per_Poverty, poor_health_percent, Index_Asian, Index_Latino
#PCA cluster 6
c6_PCA <- pca(cluster6a[2:18],transposed = T)
screeplot(c6_PCA)
findElbowPoint(c6_PCA$variance)## PC10
## 10biplot(c6_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 6')## Warning: Removed 1 rows containing missing values (geom_segment).## Warning: Removed 1 rows containing missing values (geom_label_repel).
plotloadings(c6_PCA, title="Cluster 6")## -- variables retained:## Index_Latino, Xstreamlengthimpaired, Gini, Bachelor_Over_25, Per_Avg_Land_Cov, non_migration, poor_health_percent, Index_Asian, FOOD_INS, Index_Black
#PCA cluster 7
c7_PCA <- pca(cluster7a[2:18],transposed = T)
screeplot(c7_PCA)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
findElbowPoint(c7_PCA$variance)## PC6
## 6biplot(c7_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 7')## Warning: Removed 2 rows containing missing values (geom_segment).## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c7_PCA, title="Cluster 7")## -- variables retained:## Xstreamlengthimpaired, AQI_Good, Bachelor_Over_25, Per_Avg_Land_Cov, VIO_CRIME, Index_Latino
#PCA cluster 8
c8_PCA <- pca(cluster8a[2:18],transposed = T)
screeplot(c8_PCA)
findElbowPoint(c8_PCA$variance)## PC9
## 9biplot(c8_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 8')## Warning: Removed 2 rows containing missing values (geom_segment).
## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c8_PCA, title="Cluster 8")## -- variables retained:## GHG_Percap, Index_Black, non_migration, Bachelor_Over_25, Per_Sev_Hous, UNEMPLOY, Index_Asian, Z_Water_Index, Index_Latino
#PCA cluster 9
c9_PCA <- pca(cluster9a[2:18],transposed = T)
screeplot(c9_PCA)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).
findElbowPoint(c9_PCA$variance)## PC9
## 9biplot(c9_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 9')## Warning: Removed 4 rows containing missing values (geom_segment).## Warning: Removed 4 rows containing missing values (geom_label_repel).
plotloadings(c9_PCA, title="Cluster 9")## -- variables retained:## Index_Asian, non_migration, FOOD_INS, Index_Black, Bachelor_Over_25, Gini, Per_Avg_Land_Cov, Per_Sev_Hous
#PCA cluster 10
c10_PCA <- pca(cluster10a[2:18],transposed = T)
screeplot(c10_PCA)
findElbowPoint(c10_PCA$variance)## PC9
## 9biplot(c10_PCA, showLoadings = T, showLoadingsNames = T, hline=0,vline=0, title='Cluster 10')## Warning: Removed 2 rows containing missing values (geom_segment).## Warning: Removed 2 rows containing missing values (geom_label_repel).
plotloadings(c10_PCA, title="Cluster 10")## -- variables retained:## Index_Latino, non_migration, Gini, Bachelor_Over_25, Per_Poverty, AQI_Good, Per_Sev_Hous, Xstreamlengthimpaired, Index_Black, Per_Avg_Land_Cov
Dimension reduction
#run discriminant analysis on clusters
DR <- MclustDR(dat_mc10, lambda = 1) #setting lambda to 1 gives most separating directions
summary(DR)## -----------------------------------------------------------------
## Dimension reduction for model-based clustering and classification
## -----------------------------------------------------------------
##
## Mixture model type: Mclust (VVE, 10)
##
## Clusters n
## 1 50
## 2 59
## 3 53
## 4 62
## 5 150
## 6 129
## 7 51
## 8 66
## 9 142
## 10 124
##
## Estimated basis vectors:
## Dir1 Dir2 Dir3 Dir4 Dir5
## AQI_Good -0.176001 0.0167446 0.2544551 -0.0041113 -0.199672
## Bachelor_Over_25 0.072840 -0.2897267 0.3135171 -0.5609009 0.421081
## Per_Poverty 0.022156 0.2216694 0.6106306 -0.1721047 0.485515
## Gini -0.018424 -0.1870715 -0.0774274 -0.0713467 -0.052922
## non_migration 0.122296 0.0306189 -0.1004147 0.4486234 0.108786
## Per_Sev_Hous -0.096523 -0.2903113 -0.5171334 -0.2399188 0.140611
## Xstreamlengthimpaired 0.474859 0.2044703 -0.1971267 0.1582116 0.240913
## Per_Avg_Land_Cov -0.253104 -0.3671406 -0.0276019 0.2069777 -0.338410
## poor_health_percent -0.160339 -0.3429140 -0.0287426 0.2401105 0.097547
## Z_Water_Index -0.751602 0.6152728 -0.1207784 0.0101883 0.270779
## Index_Black 0.012783 0.0364651 -0.1388193 -0.0924190 0.063046
## Index_Asian 0.030941 -0.1078755 0.1763719 0.0492666 0.086927
## Index_Latino 0.053899 0.0566457 -0.0075445 0.2151389 -0.028245
## GHG_Percap -0.134751 -0.0033918 0.1659712 0.0513501 0.460770
## UNEMPLOY 0.016303 -0.0286237 -0.2089966 0.0738562 0.156282
## FOOD_INS 0.113812 -0.2359471 0.0577083 0.3324095 -0.089733
## VIO_CRIME -0.151807 -0.0315984 -0.0236525 0.3006170 0.013413
## Dir6 Dir7 Dir8 Dir9
## AQI_Good 0.247547 0.396043 0.0874908 -0.0832171
## Bachelor_Over_25 -0.170414 -0.031725 -0.3655544 0.0382844
## Per_Poverty -0.064981 0.201735 -0.6034751 -0.4234419
## Gini -0.079564 -0.351825 0.2765681 -0.1192998
## non_migration 0.292492 -0.165084 -0.1740190 -0.2571548
## Per_Sev_Hous 0.351465 0.364975 0.3325524 -0.2378273
## Xstreamlengthimpaired -0.438677 0.568038 -0.0946994 -0.0208607
## Per_Avg_Land_Cov -0.188715 0.293543 -0.4090132 0.0051191
## poor_health_percent -0.296444 0.183355 0.1365397 0.1702002
## Z_Water_Index -0.218556 0.071558 -0.0198185 0.0456237
## Index_Black 0.105431 -0.086784 -0.1339241 0.4102571
## Index_Asian -0.134440 -0.048924 0.0340747 -0.1862869
## Index_Latino 0.026287 -0.080645 -0.1474826 -0.1506777
## GHG_Percap 0.322789 0.095990 -0.0615401 0.1927214
## UNEMPLOY 0.388493 0.056520 -0.0461231 0.2111648
## FOOD_INS -0.116806 -0.038561 0.1871896 0.5775028
## VIO_CRIME -0.166452 -0.198262 0.0084581 -0.0539742
##
## Dir1 Dir2 Dir3 Dir4 Dir5 Dir6 Dir7
## Eigenvalues 0.94731 0.67751 0.50022 0.43587 0.32617 0.18446 0.13214
## Cum. % 28.77575 49.35584 64.55074 77.79092 87.69873 93.30189 97.31587
## Dir8 Dir9
## Eigenvalues 0.06914 0.019223
## Cum. % 99.41608 100.000000#plot eigenvalues
plot(c(1:17),DR$evalues,ylab="eigenvalues", xlab="component", type="o", pch=16)
#plot all pairs of combinations
plot(DR, what='pairs',symbols=c("1","2","3","4","5","6","7","8","9","10"),
colors=c("red","blue","green","goldenrod2","violet","brown","black","grey30","slateblue2","seagreen3"))
#plot directions for variables
Directions <- data.frame(d1=DR$basis[,1],d2=DR$basis[,2],
d3=DR$basis[,3],d4=DR$basis[,4],
d5=DR$basis[,5],d6=DR$basis[,6],
d7=DR$basis[,7],d8=DR$basis[,8],
d9=DR$basis[,9],var=rownames(DR$basis))
par(mfrow=c(3,3))
dotchart(Directions$d1, labels=Directions$var, pch=16, main="Dir1", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d2, labels=Directions$var, pch=16, main="Dir2", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d3, labels=Directions$var, pch=16, main="Dir3", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d4, labels=Directions$var, pch=16, main="Dir4", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d5, labels=Directions$var, pch=16, main="Dir5", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d6, labels=Directions$var, pch=16, main="Dir6", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d7, labels=Directions$var, pch=16, main="Dir7", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d8, labels=Directions$var, pch=16, main="Dir8", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
dotchart(Directions$d9, labels=Directions$var, pch=16, main="Dir9", xlim=c(-0.6,0.6))
abline(v=0, lwd=2)
abline(v=c(0.1,-0.1),lty=2)
par(mfrow=c(1,1))